subroutine calceaab(iea, ieb, tmpe, tmpv,  i, j, k,  &
                  irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb) 
use mod_iop 
use mod_t
implicit none
integer, intent(in) :: vrtina(8), popina(8)
integer, intent(in) :: vrtinb(8), popinb(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*8, intent(inout):: iea(*), ieb(*), tmpe(*), tmpv(*)
!
if(prec==8.and.igpu==0) then 
   call calceaab_d(iea, ieb, tmpe, tmpv, i, j, k, &
               irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb)
elseif(prec==4.and.igpu==0) then 
   call calceaab_s(iea, ieb, tmpe, tmpv, i, j, k, &
               irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb)
elseif(prec==8.and.igpu==1) then 
   call calceaab_dg(iea, ieb, tmpe, tmpv, i, j, k, &
               irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb)
elseif(prec==4.and.igpu==1) then 
   call calceaab_sg(iea, ieb, tmpe, tmpv, i, j, k, &
               irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb)
endif
return
end
subroutine calceaab_d(iea, ieb, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtina(8), popina(8)
integer, intent(in) :: vrtinb(8), popinb(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*8, intent(inout):: iea(*), ieb(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab(8), ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffoa(8), ioffva(8), ioffout, idxout
integer ioffob(8), ioffvb(8)
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
do irrepab = 1, nirrep 
   inumab(irrepab) = irpdso(irrepab, vrtina, vrtina, 0) 
enddo
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffoa(1) = 0
ioffob(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffoa(irrepi2 + 1) = ioffoa(irrepi2) + popina(irrepi2)
   ioffob(irrepi2 + 1) = ioffob(irrepi2) + popinb(irrepi2)
enddo
!
ioffva(1) = ioffoa(nirrep) + popina(nirrep)
ioffvb(1) = ioffob(nirrep) + popinb(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffva(irrepi2 + 1) = ioffva(irrepi2) + vrtina(irrepi2)
   ioffvb(irrepi2 + 1) = ioffvb(irrepi2) + vrtinb(irrepi2)
enddo
!
idxi = ioffoa(irrepi) + i 
idxj = ioffoa(irrepj) + j 
idxk = ioffob(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtina, vrtina, 0)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtina(irrepa)
      numb = vrtina(irrepb)
      numc = vrtinb(irrepc)
      ioffa = ioffva(irrepa)
      ioffb = ioffva(irrepb)
      ioffc = ioffvb(irrepc)
      ioffout = isymoffso(irrepc, irrepabc, inumab, vrtinb, 1) - 1 & 
              + isymoffso(irrepb, irrepab, vrtina, vrtina, 0) - 1
      if(irrepa>irrepb) cycle
      if(irrepa==irrepb) then
!$omp parallel do
          do c = 1, numc
          do b = 2, numb 
          do a = 1, b-1
              idxout = ioffout + (c-1)*numab + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$omp end parallel do
       elseif(irrepa<irrepb) then
!$omp parallel do
          do c = 1, numc
          do b = 1, numb
          do a = 1, numa
             idxout = ioffout + (c-1)*numab+(b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$omp end parallel do
       endif
    enddo
enddo
return
end
subroutine calceaab_s(iea, ieb, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtina(8), popina(8)
integer, intent(in) :: vrtinb(8), popinb(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*4, intent(inout):: iea(*), ieb(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab(8), ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffoa(8), ioffva(8), ioffout, idxout
integer ioffob(8), ioffvb(8)
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
do irrepab = 1, nirrep 
   inumab(irrepab) = irpdso(irrepab, vrtina, vrtina, 0) 
enddo
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffoa(1) = 0
ioffob(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffoa(irrepi2 + 1) = ioffoa(irrepi2) + popina(irrepi2)
   ioffob(irrepi2 + 1) = ioffob(irrepi2) + popinb(irrepi2)
enddo
!
ioffva(1) = ioffoa(nirrep) + popina(nirrep)
ioffvb(1) = ioffob(nirrep) + popinb(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffva(irrepi2 + 1) = ioffva(irrepi2) + vrtina(irrepi2)
   ioffvb(irrepi2 + 1) = ioffvb(irrepi2) + vrtinb(irrepi2)
enddo
!
idxi = ioffoa(irrepi) + i 
idxj = ioffoa(irrepj) + j 
idxk = ioffob(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtina, vrtina, 0)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtina(irrepa)
      numb = vrtina(irrepb)
      numc = vrtinb(irrepc)
      ioffa = ioffva(irrepa)
      ioffb = ioffva(irrepb)
      ioffc = ioffvb(irrepc)
      ioffout = isymoffso(irrepc, irrepabc, inumab, vrtinb, 1) - 1 & 
              + isymoffso(irrepb, irrepab, vrtina, vrtina, 0) - 1
      if(irrepa>irrepb) cycle
      if(irrepa==irrepb) then
!$omp parallel do
          do c = 1, numc
          do b = 2, numb 
          do a = 1, b-1
              idxout = ioffout + (c-1)*numab + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$omp end parallel do
       elseif(irrepa<irrepb) then
!$omp parallel do
          do c = 1, numc
          do b = 1, numb
          do a = 1, numa
             idxout = ioffout + (c-1)*numab+(b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$omp end parallel do
       endif
    enddo
enddo
return
end
subroutine calceaab_dg(iea, ieb, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtina(8), popina(8)
integer, intent(in) :: vrtinb(8), popinb(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*8, intent(inout),device:: iea(*), ieb(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab(8), ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffoa(8), ioffva(8), ioffout, idxout
integer ioffob(8), ioffvb(8)
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
do irrepab = 1, nirrep 
   inumab(irrepab) = irpdso(irrepab, vrtina, vrtina, 0) 
enddo
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffoa(1) = 0
ioffob(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffoa(irrepi2 + 1) = ioffoa(irrepi2) + popina(irrepi2)
   ioffob(irrepi2 + 1) = ioffob(irrepi2) + popinb(irrepi2)
enddo
!
ioffva(1) = ioffoa(nirrep) + popina(nirrep)
ioffvb(1) = ioffob(nirrep) + popinb(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffva(irrepi2 + 1) = ioffva(irrepi2) + vrtina(irrepi2)
   ioffvb(irrepi2 + 1) = ioffvb(irrepi2) + vrtinb(irrepi2)
enddo
!
idxi = ioffoa(irrepi) + i 
idxj = ioffoa(irrepj) + j 
idxk = ioffob(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtina, vrtina, 0)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtina(irrepa)
      numb = vrtina(irrepb)
      numc = vrtinb(irrepc)
      ioffa = ioffva(irrepa)
      ioffb = ioffva(irrepb)
      ioffc = ioffvb(irrepc)
      ioffout = isymoffso(irrepc, irrepabc, inumab, vrtinb, 1) - 1 & 
              + isymoffso(irrepb, irrepab, vrtina, vrtina, 0) - 1
      if(irrepa>irrepb) cycle
      if(irrepa==irrepb) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 2, numb 
!$acc loop independent
          do a = 1, b-1
              idxout = ioffout + (c-1)*numab + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$acc end kernels 
       elseif(irrepa<irrepb) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 1, numb
!$acc loop independent
          do a = 1, numa
             idxout = ioffout + (c-1)*numab+(b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$acc end kernels 
       endif
    enddo
enddo
return
end
subroutine calceaab_sg(iea, ieb, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtina, vrtinb, popina, popinb) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtina(8), popina(8)
integer, intent(in) :: vrtinb(8), popinb(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*4, intent(inout),device:: iea(*), ieb(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab(8), ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffoa(8), ioffva(8), ioffout, idxout
integer ioffob(8), ioffvb(8)
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
do irrepab = 1, nirrep 
   inumab(irrepab) = irpdso(irrepab, vrtina, vrtina, 0) 
enddo
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffoa(1) = 0
ioffob(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffoa(irrepi2 + 1) = ioffoa(irrepi2) + popina(irrepi2)
   ioffob(irrepi2 + 1) = ioffob(irrepi2) + popinb(irrepi2)
enddo
!
ioffva(1) = ioffoa(nirrep) + popina(nirrep)
ioffvb(1) = ioffob(nirrep) + popinb(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffva(irrepi2 + 1) = ioffva(irrepi2) + vrtina(irrepi2)
   ioffvb(irrepi2 + 1) = ioffvb(irrepi2) + vrtinb(irrepi2)
enddo
!
idxi = ioffoa(irrepi) + i 
idxj = ioffoa(irrepj) + j 
idxk = ioffob(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtina, vrtina, 0)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtina(irrepa)
      numb = vrtina(irrepb)
      numc = vrtinb(irrepc)
      ioffa = ioffva(irrepa)
      ioffb = ioffva(irrepb)
      ioffc = ioffvb(irrepc)
      ioffout = isymoffso(irrepc, irrepabc, inumab, vrtinb, 1) - 1 & 
              + isymoffso(irrepb, irrepab, vrtina, vrtina, 0) - 1
      if(irrepa>irrepb) cycle
      if(irrepa==irrepb) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 2, numb 
!$acc loop independent
          do a = 1, b-1
              idxout = ioffout + (c-1)*numab + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$acc end kernels 
       elseif(irrepa<irrepb) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 1, numb
!$acc loop independent
          do a = 1, numa
             idxout = ioffout + (c-1)*numab+(b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(iea(idxi) + iea(idxj) + ieb(idxk) &
                          -iea(idxa) - iea(idxb) - ieb(idxc))
!              write(6,'(A,I4, I4, I4, I4,I4)')'idxout', irrepi, irrepj, irrepk, ioffout, idxout
          enddo
          enddo
          enddo
!$acc end kernels 
       endif
    enddo
enddo
return
end
